Molecular Ecology Resources

Abstract Background Many ecologically important small invertebrates are preserved in ethanol, which compromises DNA integrity and severely limits their potential to be utilized in further genomic research. The lack of soil arthropod representation ultimately results in limited availability of reference genomes and population genomic inference from pooled sequencing (Pool-seq). Results We introduce PoLoCo (Pooled Low-Coverage), a reproducible workflow integrating pooled specimen draft genome assembly from ethanol-preserved specimens with population-level allele-frequency analysis. We assembled a 411 Mb draft genome of Entomobrya nivalis from pooled short reads, achieving 97.7% BUSCO completeness despite notable fragmentation. Comparison with the published reference genome showed partial genomic correspondence while also revealing reference-dependent differences in the retained SNP datasets. The application of PoLoCo to 82 retained pooled populations produced a final distance-thinned allele-frequency dataset suitable for population-scale genomic inference. Case-study execution documented the computational requirements of the major workflow components on standard HPC infrastructure, with the workflow, documentation, and software environments openly available in Zenodo. Conclusions PoLoCo demonstrates that draft genome assemblies and reproducible Pool-seq allele-frequency analyses can be performed with long-preserved sample material from the field. PoLoCo integrates reference construction, sync-based allele-frequency analysis, and benchmarking into a free and open workflow. It establishes a reproducibility benchmark for biodiversity genomics and yields a transferable workflow for ecological and evolutionary initiatives where high-quality DNA or long-read sequencing is not feasible.

Environmental DNA (eDNA) has emerged as a transformative tool for biodiversity monitoring, offering a non-invasive means to detect and identify species across various ecosystems. However, despite its potential, current eDNA methodologies have significant limitations, including limited sensitivity, restriction to a specific taxonomic group, lengthy processing times, the need for specialized laboratory conditions and equipment, as well as computational bottlenecks in equipment, connectivity, and algorithms. Overcoming these limitations is crucial for the real-time monitoring of a high number of taxa in field-based scenarios. Here we present an integrated eDNA workflow from extraction to data analysis that mitigates several of these challenges through several technical innovations. Our approach adapts the Tagsteady protocol to Nanopore sequencing, integrates multiplex PCR for simultaneous plant and animal detection, and introduces RTBioScan, a real-time analysis pipeline employing hierarchical basecalling and OTU-level consensus generation to achieve high-accuracy taxonomic assignments during sequencing. Using a single set of universal primers, we integrate metabarcoding and barcoding into a cohesive protocol, utilizing universal markers for both animals and plants (cytochrome oxidase 1 region, COI, and Internal Transcribed Spacer 2, ITS2). This comprehensive approach can detect and identify a broad spectrum of species, from individual organisms to entire communities in a single workflow.

Structural variants (SVs) can strongly influence evolutionary processes by suppressing recombination to create jointly inherited regions and by altering gene expression and the 3D organization of DNA. Recent improvements in sequencing technologies have led to an influx of high-resolution genomic data, bolstering research on SVs. When detecting and characterizing SVs, researchers face trade-offs in data and methods concerning scale, resolution, and computational intensity. Here, we evaluate such trade-offs by quantifying how data type, analytical approach, and SV attributes affect the ability to detect and characterize SVs in Timema cristinae stick insects, with a focus on inversions. We specifically evaluate inversions characterized by (i) a population genomic approach based on patterns of local population structure inferred from genotyping-by-sequencing data, (ii) pairwise alignments of eight de novo genome assemblies and SV calling with SyRI, and (iii) a pangenome constructed from the alignment of the eight genome assemblies. Our results suggest that while broad genomic regions harboring SVs can be detected by any of the approaches, the approaches differ in how they characterize inversions, especially small or complex ones. The local population structure approach is exploratory and generally only detects the largest inversions. Pairwise comparative alignments identify the greatest number of inversions, but downstream analyses are influenced by identification of shared inversions across genome pairs. Pangenome approaches are scalable and produce network graphs that describe complex SVs but require a conceptual shift in visualization and interpretation. The collective results highlight limitations and benefits of different approaches in a burgeoning research field.

During biogeographical ancestry (BGA) analysis, an unknown individual for whom ancestry informative markers (AIMs) were analyzed is often assigned to one of a few reference populations. For some individuals, such an assignment is not meaningful: the two parents of the individual may come from different populations, or the individual may come from a population not included in the reference database. In this study, we examined classification methods in combination with 104 AIMs from the VISAGE Enhanced Tool to provide one of the following three possible outcomes: (i) in-reference (IR), with an assignment to one of the reference populations; (ii) admixed (AD), with parents from two in-reference populations; and (iii) out-of-reference (OR), where the individual's ancestry is not sufficiently represented by the reference populations. We therefore compared seven classifiers in combination with class-specific acceptance thresholds for the OR decision. The best-performing method according to our selection criterion was a naive Bayes classifier over pairs of populations, which we call the recent naive Bayes (RNB) classifier. Populations held out of the reference panel were correctly rejected in $84.9\%$ of cases. We further tested the method on $800$ individuals from seven additional groups not included in the reference panel. Only $29.4\%$ of individuals from four non-admixed populations were rejected, resembling joint ancestry with reference populations. In addition, we investigated the prediction accuracy in a newly collected dataset of 28 individuals with self-reported ancestry by analyzing AIMs using massively parallel sequencing.

Polyandry, multiple mating by females, both benefits females by increases genetic variation within families and provides an opportunity for sexual selection through sperm competition and cryptic female choice. Estimating polyandry in species with large clutches is challenging because genotyping all progeny individually can be impractical. We present a likelihood-based framework for estimating the number of sires and their relative contributions from the maternal genotype and pooled sequencing of entire progeny pools. The method evaluates the probability of observed pooled offspring read counts under alternative paternal genotype mixtures, accounting for Mendelian segregation and sampling variation from finite sequencing depth. We implement this framework for both biallelic SNPs and short-read mini-haplotype markers constructed from closely linked SNPs, which provide multi-allelic information comparable to traditional microsatellite markers. Simulations show that this approach accurately recovers sire number and sire contributions across a range of conditions. We then applied the method to 275 wild-caught Drosophila melanogaster females and their offspring from a Kansas population using 3,053 biallelic SNPs and 239 four-allele mini-haplotype markers. Estimates from the two marker sets were highly concordant, and combined analyses indicate that 57% of families were sired by a single male. These results demonstrate that pooled progeny sequencing provides an efficient genomic approach for estimating polyandry in species with large families.

Global arthropod decline demands effective biodiversity monitoring strategies. However, most current monitoring approaches do not provide an exhaustive picture of arthropod community structure. In particular, biotic interactions and temporal patterns of biodiversity change are still poorly understood. Here we explore the possibility of addressing these two shortfalls using spiders, one of the most important predators of terrestrial arthropods, as natural samplers for arthropod community DNA. We conducted several experiments comparing the recovered community composition between spider gut contents and traditional monitoring methods. Additionally, we used archived spiders from long-term collections to assess the preservation of prey DNA over time. Spiders proved to be highly efficient natural DNA samplers, with gut content metabarcoding revealing similar community composition and alpha- and beta-diversity like traditional methods. Unique arthropod taxa were detected by spider gut contents and traditional methods respectively, indicating that these two approaches can serve as valuable complements to each other. Well-preserved archived spiders effectively reconstructed historical diets up to 20 years back in time, showing that historical spider collections can act as time capsules of the wider diversity that occurred in the spider's habitat. Spiders as natural samplers overcome critical shortfalls in biodiversity monitoring and contribute to our future understanding of community assembly across space and time.

At large scales, our ability to compare genetic diversity across species is often challenged by the difficulty of extracting range-wide sequence data from online databases. Using mitochondrial DNA from amphibians and reptiles spanning across six regions: Brazilian Atlantic Forest, Iberian Peninsula, northwestern Africa, Arabian Peninsula, Madagascar and Australian Monsoonal Tropics, we develop and implement a framework to overcome common limitations in macrogenetics. We compiled and manually curated a spatially explicit genetic dataset from public repositories, quantified the spatial representativeness of genetic sampling across species' ranges and developed a novel interpolation framework to estimate and map species-specific patterns of intraspecific genetic diversity (ISD) while accounting for heterogeneous spatial sampling. We estimated mean nucleotide diversity across species to produce regional maps for each genetic marker. Our final dataset comprised 182 amphibian and 239 reptile species, encompassing 184 genera from 54 families. Range representativeness of genetic sampling was generally high, despite some variation across study regions. Spatial patterns of ISD revealed geographic heterogeneity, and were consistent across the best-sampled mitochondrial markers. By combining automated and manual georeferencing, quantifying range representativeness, and accounting for bias in spatial interpolations, our pipeline is marker-agnostic and directly applicable to nuclear and genomic data as they become available.

Detecting dugongs at low densities remains a major challenge for conservation monitoring, particularly in remote or turbid coastal habitats where conventional surveys are difficult. Environmental DNA (eDNA) offers a complementary approach, but robust species-specific tools are needed before it can be confidently applied to dugong surveillance. We developed and analytically validated a quantitative PCR TaqMan assay for detecting dugong (Dugong dugon) eDNA and evaluated its performance across controlled aquarium conditions in Australia and Abu Dhabi, and wild populations in northern Australia and the Arabian/Persian Gulf. The assay detected dugong eDNA in aquarium and field samples, with Sanger sequencing confirming amplification of dugong mitochondrial DNA from positive field detections. Our assay provides a robust and geographically transferable tool for targeted dugong eDNA detection, performing reliably across contrasting coastal environments and widely separated populations. This approach can strengthen local-scale surveillance of low-density or difficult-to-detect dugong populations.

Background Comparisons of single-cell and single-nucleus RNA sequencing (scRNA-seq and snRNA-seq) data have been described in some mammalian tissues and, subsequently, in Drosophila, but remain unexplored in most invertebrate lineages. The xenacoelomorphs occupy key phylogenetic positions, yet they differ anatomically from mammals. They have a reduced extracellular matrix, high-salt body fluid, and no circulatory system. Despite these differences, they possess a well-developed nervous system. One of the xenacoelomorphs, the photosymbiotic acoel (Praesagittifera naikaiensis) also harbours symbiotic Tetraselmis algae, whose RNA can be co-captured with host RNA. Results We compared scRNA-seq and snRNA-seq data from whole P. naikaiensis specimens. Both methods yielded high-quality data with comparable gene detection but a larger share of scRNA-seq reads derived from symbiotic algae. Gene-level analyses revealed that neural genes were enriched in snRNA-seq relative to non-neural genes. Cross-method label transfer and integration-based validation identified six snRNA-seq clusters, including some neural populations, that lacked a clear scRNA-seq counterpart. In contrast, three cell populations, including muscle and metabolically active clusters, were reciprocally validated as captured by both methods. Conclusions Our results show that key snRNA-seq advantages, particularly the enhanced recovery of neural transcripts, are recapitulated in our dataset, which is consistent with previous reports in mammals. Furthermore, snRNA-seq reduces symbiont-derived reads and recovers several cell populations underrepresented in scRNA-seq. These findings provide practical guidance for cell atlas construction in non-model, symbiotic invertebrates.

A complete reference genome for species with chromosomally-determined separate sexes should contain scaffolds for all sex chromosome homologs. However, sex chromosomes present distinct computational challenges compared to autosomes. Here we present a k-mer based analysis that utilizes whole-genome sequencing of a few sex-identified isolates: Cytogenetics-By-Sequencing (CBS) tools. Unlike other approaches that typically address one aspect of the sex chromosomes, CBS-tools strives to guide users from the discovery of the heterogametic sex through identifying the sex-determination region (SDR). The core of CBS-tools is automated quantification of sex-specific k-mers in order to predict the heterogametic sex. Using publicly-available datasets, CBS-tools correctly identified the known sex-system of the 31 species tested. Additionally, we used these k-mers to verify and correct phasing of sex chromosomes between haplotypes in species representing different sex-systems. Finally, we used these k-mers to delimit the SDR boundary using an interactive web platform. CBS-tools was developed with previously unexplored sex chromosome systems in mind, but is also suitable for well-examined sex chromosome pairs.

Understanding the demographic processes shaping population connectivity is essential for conservation. Yet disentangling the effects of dispersal and local effective density remains challenging because traditional genetic approaches integrate signals over evolutionary timescales. Identity-by-descent (IBD) segments provide opportunities to reconstruct recent demographic connectivity over timescales relevant to contemporary conservation challenges, but their reliability under realistic ecological conditions remains poorly evaluated. Here, we assess the ability of two promising spatiall IBD-based methods, IBD-Analysis and MAPS, to estimate effective density and dispersal variance. Using a simulation-inference framework combining spatially explicit simulations, IBD segment extraction, and demographic inference, we evaluated both methods across a wide range of demographic, ecological and methodological scenarios. Both approaches reliably recovered demographic parameters under homogeneous landscapes, with IBD-Analysis providing accurate estimates of global effective density and dispersal across a broad range of conditions. Spatial heterogeneity emerged as the main driver of inference performance. When demographic parameters varied across space, IBD-Analysis integrated local variation into biased global estimates, whereas MAPS captured spatial patterns but showed substantial variability and generated spurious spatial structure in homogeneous landscapes. Methodological factors had comparatively weaker effects than ecological complexity, highlighting the importance of accounting for landscape structure when applying IBD-based demographic inference. These results highlight the potential of IBD-based approaches to move beyond descriptive patterns of genetic connectivity toward mechanistic inference of demographic processes. By distinguishing the contributions of dispersal and density, these methods can improve conservation decisions by linking genomic patterns to key demographic parameters relevant for conservation.

Generating high-quality genome assemblies for small animals with large genomes is complex due to their small body size, DNA contamination, and repetitive elements. Ticks exemplify these complexities, while also being a global health threat to humans, domestic animals, and wildlife. Advances in long-read sequencing platforms now make it feasible to obtain large amounts of raw sequence data from individual specimens, but challenges remain. Key genome workflow challenges include error correction, assembly, transposable element annotation, and chromosome assignment and synteny. Here we examine three individual Dermacentor reticulatus ticks using deep Oxford Nanopore sequencing. Comparing and contrasting bioinformatic tools for raw read and assembly manipulation allows us to identify the parameters that provide a high-quality haploid genome assembly among 38 assemblies. We find Dorado error corrected raw read data from approximately two flow cells are needed, but limiting assembly coverage to 40x produces the highest quality assemblies. Then we examined the effects of optimal workflow parameters on downstream analyses of gene synteny and manual transposable element annotation of the three tick assemblies in comparison to an independent assembly of a fourth D. reticulatus individual from a different country. Gene synteny analysis allows chromosome assignment of scaffolds and transposable element identification was improved markedly with manual curation. Finally, we compared the genetic variation of D. reticulatus across two populations and found similar genetic diversity. Our study provides a clear workflow to obtain high-quality assemblies from Oxford Nanopore sequences and genetic characterization of a species with a large, complex, and repetitive genome.

Background: Amplicon sequencing of the 16S rRNA gene is widely used to characterise microbial communities, but the performance of commonly used clustering and denoising pipelines has not been systematically evaluated for highly diverse environmental datasets. We therefore benchmarked four established clustering and denoising pipelines, VSEARCH cluster_size, UNOISE as implemented in VSEARCH, Swarm, and DADA2, using simulated microbial communities with known ground-truth compositions and evaluated how performance varied with species richness, abundance unevenness, sequencing depth, and minimum abundance threshold. Results: Differences among pipelines were small at low species richness but became greater as species richness increased. High average cluster purity did not necessarily correspond to accurate species-level recovery, as species could be split across multiple clusters or recovered incompletely. Among the pipelines, UNOISE may be preferable when high species representation and cluster purity are prioritised, while also showing strong reconstruction of relative abundance profiles, albeit with extensive species splitting and many unclustered reads. DADA2 showed similarly strong reconstruction of relative abundance profiles and less species splitting but represented fewer species and had lower cluster purity at higher richness. Swarm provided a balanced compromise between species representation and splitting, whereas cluster_size may be useful when limiting species splitting is a priority, despite weaker reconstruction of relative abundance profiles. Increasing the minimum abundance threshold reduced species splitting but also reduced the number of clusters and perfect clusters and, at higher thresholds, decreased concordance with ground-truth relative abundance profiles. These effects were more pronounced at lower sequencing depths. Analysis of seafloor sediment samples also showed threshold-dependent sequence loss, including the loss of sequences repeatedly detected across replicate samples. Conclusions: Pipeline performance in highly diverse 16S rRNA amplicon datasets varies with dataset characteristics, evaluation criteria, and parameter settings. Pipeline selection should therefore reflect dataset characteristics and the analytical objective rather than rely on a single measure of performance, and minimum abundance thresholds should be evaluated in relation to sequencing depth rather than applied as fixed defaults. These findings emphasise the importance of reconsidering established analytical practices and transparently reporting bioinformatic settings as 16S rRNA amplicon sequencing is increasingly applied to highly diverse environmental microbial communities.

The hybrid grouper 'Yushuban' (Cromileptes altivelis [female] x Epinephelus fuscoguttatus [male]) has shown promising performance in aquaculture trials, yet no complete reference genome has been available for this hybrid. We generated two gapless, telomere-to-telomere (T2T) haplotype-resolved genome assemblies of a single individual by integrating short-read, PacBio HiFi, Oxford Nanopore ultralong-read and Hi-C data. The maternal (Yushu_Calti_1.0) and paternal (Yushu_Efusc_1.0) haplotypes are 1,077.33 Mb and 1,083.36 Mb in length, with N50 values of 45.46 Mb and 45.55 Mb, respectively. Across all platforms, over 99.4% of reads mapped to each assembly, and Merqury-based quality values reached 57.7413 and 59.2068 (PacBio). BUSCO recovered 99.78% and 99.75% of complete core genes. We predicted 27,248 and 27,566 protein-coding genes, of which 98.52% and 99.23% were functionally annotated, and annotated telomeric motifs and repetitive elements (46.12% and 47.18% of each genome). These assemblies fill a gap in grouper genomic resources and provide a foundation for genetic dissection, parental genome comparison and molecular breeding of this hybrid.

Understanding floral visitor dynamics along elevational gradients is essential for assessing potential implications for pollen limitation and population genetic processes in fragmented populations of cliff-dwelling alpine endemics. Physoplexis comosa (Campanulaceae) is a chasmophytic paleoendemic restricted to vertical rock faces on dolomitic and limestone bedrock from 200 to 2,800 m a.s.l. The species exhibits a specific floral morphology with fused corolla apices, basal openings for nectar access, and secondary pollen presentation, yet its functional floral visitor assemblage has remained poorly characterized across the species' elevational range. In this pilot study, we tested a two-tier camera-monitoring setup combining action cameras for diurnal monitoring with Raspberry Pi-based systems equipped with IR illumination for combined diurnal and nocturnal recording. We conducted monitoring across five elevational sites (1,480 to 2,400 m a.s.l.) in the Eastern Alps. Recorded video material was downsampled and annotated using the Addax platform alongside manual cross-validation to record physical contact events. Preliminary results confirm the technical feasibility of continuous monitoring in steep terrain and suggest possible elevational variation in diurnal visitors, alongside nocturnal activity by Geometridae moths. Camera-based observations provided higher functional resolution and revealed a different visitor assemblage than previous eDNA metabarcoding trials. This proof-of-concept establishes a methodological foundation for forthcoming multi-day monitoring campaigns and complements ongoing studies of the mating system and population-genetic analyses.

Motivation: Optimum contribution selection (OCS) maximizes genetic gain under a constraint on relatedness. Genomic relationship matrices are dense, so the quadratic term costs O(n^3) operations, and existing OCS software assumes pedigree sparsity. Results: We replace the relationship matrix by its leading eigenpairs, computed by randomized SVD, and solve OCS in factorized form. On Norway spruce population (n=5,525), ranks 10--30 kept the exact objective within 1% of the optimum at 110--330-fold speedups; on the QTL-MAS 2010 (n=3,226), ranks up to 100 were needed (24-211 fold). The required rank rises as the diversity constraint loosens. Gain alone is a misleading accuracy check, since coancestry exceeded the optimum by up to 28%, so we choose the rank on the exact objective, which can be checked without a dense solve. The speedup comes from computing only the leading eigenpairs. We extend the factorization to robust OCS (18-fold). Availability and implementation: Julia implementation, and examples are available at https://github.com/jonhar97/fastocs under the MIT license (archived at Zenodo, DOI 10.5281/zenodo.23214804).

Shotgun sequencing of DNA mixtures has the potential to become an important tool in forensic genetics, particularly for degraded samples or complex mixtures. However, the analysis of shotgun-sequenced mixtures is challenged by several confounding factors, including sequencing errors, population structure, the potential presence of relatives, linkage, and linkage disequilibrium. Additionally, the statistical methods currently available for such data are limited. Here, we develop likelihood ratios that summarize the evidence for the presence of a suspects DNA in a mixture, assuming genotypes of the suspect and the victim are available, while accounting for DNA from unknown contributors modeled either as general population contamination or as a single contaminating individual. We show that these likelihood ratios are not robust to model assumptions and can falsely support the presence of the suspects DNA when the true contributor is even a distant relative. We therefore propose a combined evidence score, {Lambda}: the likelihood ratio minimized over contamination models, reported only when the evidence favors the suspect over a hypothesized sibling by at least a factor of C (default C = 10) under at least one contamination model, and otherwise set to 1 (log {Lambda} = 0). In simulations based on 1000 Genomes data, {Lambda} retains nearly all power to detect true contributors at mixture proportions as low as 5% while eliminating the false positives produced by model mis-specification, population structure, and by untested relatives. The methods are implemented in the open-source program DNAMIXTURE, providing the forensic and population genetics community a resource for robustly identifying small individual contributions to a DNA mixture.

Freshwater biodiversity is declining at alarming rates globally, yet our capacity to monitor it depends heavily on the methods used for biodiversity assessment. Stream benthic macroinvertebrates are widely used as bioindicators, but traditional morphological identification is labor-intensive and limited by taxonomic expertise, while DNA metabarcoding offers a promising but incompletely evaluated alternative. Here, we compared biodiversity patterns detected by four macroinvertebrate sampling and identification methods (morphological identification of benthic samples, morphological identification of combined benthic and stream-edge samples, eDNA metabarcoding of bulk benthic samples, and eDNA metabarcoding of streamwater) across 62 sites in seven watersheds spanning the southern continental United States. We compared methods across three facets of biodiversity (gamma diversity, alpha diversity, and taxonomic composition) and three taxonomic levels (family, genus, and species). We found that which method detected more taxa shifted with taxonomic resolution: morphological methods detected more families and genera, while metabarcoding methods detected more species. Bulk benthic metabarcoding detected more families and genera than water eDNA metabarcoding, consistent with the "watered-down biodiversity" effect of DNA dilution in water samples, but the two metabarcoding methods did not differ significantly in species richness. All four methods detected largely distinct suites of taxa, with metabarcoding methods yielding more unique taxa than morphological methods. Notably, taxa uniquely detected by water eDNA were dominated by soft-bodied organisms (oligochaetes, leeches, earthworms) and water-column-associated taxa, while benthic eDNA uniquely detected hard-bodied EPT insects. We attribute this pattern to fundamental differences in eDNA shedding, transport, and deposition dynamics across organism types. These results demonstrate that no single method captures a complete picture of stream macroinvertebrate biodiversity, and that method choice should be guided by the taxonomic resolution and community components most relevant to the goals of the bioassessment program.

Automated monitoring camera systems offer an efficient approach for quantifying pollinator biodiversity and plant-pollinator interactions, but financial and logistical constraints limit the number of flowering plant species that can be monitored. Using a large European database of plant-pollinator networks, we evaluated whether plant subsampling can capture key metrics of interest (i.e. pollinator richness and network structure). We compared abundance-based, flower-shape-informed, phylogenetically informed and random plant sampling. Abundance-based sampling consistently outperformed random sampling, while adding flower shape provided little additional benefit and phylogenetic selection performed similarly to random sampling. Overall, monitoring 12-15 flowering species was sufficient to characterize key network properties across varying community sizes. Performance of abundance-based plant sampling declined in species-rich communities and when rare but highly attractive plants were present in the community, while a specific subset of pollinator species was consistently missed even by the best sampling strategy. Based on their relative contribution to interactions within each network, only 12.6% of pollinator species accounted for 95% of interactions, suggesting that AI classifiers could prioritize a relatively small subset of species. Our results provide practical guidance for designing efficient camera-based pollinator monitoring schemes at large spatial scales.

Species delimitation increasingly relies on multiple sources of evidence, but most frameworks are still limited to genetics and morphology or analyze data types separately and compare results qualitatively. We present delimSOM 2.0, an R package that leverages unsupervised machine learning through multilayer self-organizing maps. Any kind of input data, including genomic, phenotypic, ecological, and spatial, is analyzed jointly, retaining each data type as a separate layer and balancing their contributions through layer-specific weighting. Our package provides functions for preprocessing, replicate training, multivariate clustering, hyperparameter optimization, model-quality assessment, and variable- and layer-importance analyses. We evaluated our approach using two simulation sets and eight empirical studies spanning diverse taxa and data types. Overall, performance was strongest with larger sample sizes (>10 individuals per lineage) and modest numbers of lineages (K < 5), but unlike most species delimitation methods, it can reliably infer K = 1. Across two-population demographic simulations of genetic data, our method performed similarly to the well-established DAPC across gene-flow and sampling regimes, was more conservative than STRUCTURE and sNMF in inferring K = 2 and showed mixed support for K = 1 and K = 2 at intermediate differentiation levels and with gene flow. Crucially, our method is robust to missing data and variation in hyperparameters. Empirical examples showed that genetic data were typically most influential, but phenotypic, ecological, and spatial layers often contributed substantial complementary information. By combining diverse data types into a single unsupervised analysis, delimSOM 2.0 enables fully integrative species delimitation, reduces over-splitting, and advances efforts toward more objective and standardized taxonomy and speciation biology.